Double Pendulum.

01A pendulum on a pendulum

Let a pen swing loosely between two fingers: that is a pendulum, and since Galileo we have known how to predict it. Now hang a second pen from the tip of the first.

You carry a double pendulum with you every day. Let your arm hang and swing it from the shoulder: the upper arm swings about the shoulder, the forearm about the elbow, and the hand follows a path that neither joint alone would draw. A gymnast turning on the high bar, a nunchaku, a crane hook hanging from its cable are all the same machine: two masses, two rods, one gravitational field.

Released from a small angle, the double pendulum rocks gently and regularly, the two arms passing energy back and forth like two children on linked swings. For the high-angle start shown here, the lower arm whips over the top, stops and reverses in an irregular pattern. Nothing random is involved. The equations we are about to derive fix the entire future from the starting state, and the computer solves them so accurately that over the first twenty seconds the energy stays fixed to better than one part in ten million. Yet two starts that differ by a thousandth of a radian, one millimetre at the end of a one-metre rod, end up in unrelated places within five seconds.

This pairing of exact laws with practically unpredictable outcomes is called deterministic chaos. It is the same mechanism that limits weather forecasts to about two weeks: the atmosphere obeys known equations, but we can never measure today’s state precisely enough.

02The model

Physics begins by deciding what to leave out. We keep two point masses on rigid rods that weigh nothing, swinging in one vertical plane, with no friction and no air. Put the pivot at the origin with \(y\) pointing up. Reading the positions off Figure 1 with a little trigonometry,

Two angles fix the shape completely, so the system has two degrees of freedom. To know where it goes next we also need how fast each angle is changing, \(\omega_i = \mathrm{d}\theta_i/\mathrm{d}t\). The four numbers \((\theta_1, \theta_2, \omega_1, \omega_2)\) are the state: give them to the equations and the whole future follows. On the page the upper mass and rod are fixed at \(m_1 = 1\ \mathrm{kg}\), \(l_1 = 1\ \mathrm{m}\), with \(g = 9.81\ \mathrm{m/s^2}\); the sliders set the ratios \(m_2/m_1\) and \(l_2/l_1\).

03Energy

Instead of chasing forces, we will describe the pendulum by its energy. Differentiate the positions (1) with respect to time, remembering that \(\theta_i\) changes and \(l_i\) does not:

The kinetic energy is \(\tfrac12 m v^2\) for each mass. Squaring and adding the components, the mixed products combine through \(\cos\theta_1\cos\theta_2 + \sin\theta_1\sin\theta_2 = \cos(\theta_1-\theta_2)\), and the potential energy is \(mgy\) for each mass:

The last term in (3) carries the whole story. The lower mass rides on the upper one, the way a passenger walking down the aisle rides on a train. When the two rods point the same way and turn the same way, \(\cos(\theta_1-\theta_2) = 1\) and the speeds simply add, \(v_2 = l_1\omega_1 + l_2\omega_2\). When the rods are at right angles, the cosine vanishes and the two speeds add like the sides of a right triangle. Because the rods' directions keep changing, the way the two motions combine keeps changing too: each rod is forever pushing and pulling on the other.

04From energy to motion★

With Newton’s laws we would need the tension in each rod, two unknown forces that change direction and size at every instant. Lagrange found a way around them. Form the Lagrangian \(L = T - V\). Among all the ways the angles could get from a given start to a given end, the real motion is the one that makes the action \(S = \int L\,\mathrm{d}t\) stationary (Hamilton’s principle). The rod forces never appear, because they act along the rods and do no work on any motion the rods allow.

4.1The Euler–Lagrange equation

Change one angle by a small amount \(\delta\theta_i(t)\) that vanishes at the start and at the end. The action changes by \(\delta S = \int \left[ \partial_{\theta_i}L\,\delta\theta_i + \partial_{\omega_i}L\,\delta\omega_i \right]\mathrm{d}t\), with \(\delta\omega_i = \mathrm{d}(\delta\theta_i)/\mathrm{d}t\). Integrating the second term by parts moves the time derivative across; the boundary term vanishes because \(\delta\theta_i\) does:

The quantity \(p_i = \partial L/\partial\omega_i\) is the generalised momentum belonging to \(\theta_i\). For the upper angle it is \(p_1 = (m_1+m_2)\,l_1^2\,\omega_1 + m_2 l_1 l_2\,\omega_2\cos\delta\), with \(\delta = \theta_1 - \theta_2\): an angular momentum that already contains the lower mass, not simply \(m_1 l_1^2\omega_1\).

4.2The upper angle

Two kinds of derivative meet in (5). A partial derivative such as \(\partial L/\partial\omega_1\) treats \(L\) as a function of four independent numbers and varies one of them alone. The time derivative then follows the actual motion, in which all four change together; by the chain rule \(\tfrac{\mathrm{d}}{\mathrm{d}t}\cos\delta = -\sin\delta\,(\omega_1-\omega_2)\). Writing \(\alpha_i = \mathrm{d}\omega_i/\mathrm{d}t\) for the angular accelerations,

In the difference the two \(\omega_1\omega_2\) terms cancel. Dividing by \(l_1\):

4.3The lower angle

The same three steps for \(\theta_2\), with \(\partial L/\partial\omega_2 = m_2 l_2^2\,\omega_2 + m_2 l_1 l_2\,\omega_1\cos\delta\), again cancel the \(\omega_1\omega_2\) terms; dividing by \(m_2 l_2\):

4.4Solving for the accelerations

Each equation contains both accelerations, because each rod pushes on the other. A computer stepping forward in time needs each acceleration on its own, so we solve the linear system

The determinant can never vanish, so there is always exactly one answer. Cramer’s rule, then \(2\sin^2\delta = 1-\cos 2\delta\), \(2\sin\theta_2\cos\delta = \sin\theta_1 - \sin(\theta_1 - 2\theta_2)\) for \(\alpha_1\), and \(\sin\theta_1\cos\delta - \sin\theta_2 = \cos\theta_1\sin\delta\) for \(\alpha_2\) (expand \(\sin\theta_2 = \sin(\theta_1 - \delta)\); this is where the common factor \(\sin\delta\) of \(\alpha_2\) comes from), bring it to the form the page computes:

Checked against a direct numerical solve of (8) at 10 000 random states and parameter values: worst relative difference 8 × 10⁻¹⁴, which is rounding error.

05Small swings: two normal modes

Before facing the full equations, look at the gentle case. For small angles, \(\sin\theta\approx\theta\), \(\cos\delta\approx 1\), and the \(\omega^2\) terms are too small to matter. With equal masses and equal lengths \(l\), equations (6) and (7) become linear:

Linear equations like these have special solutions in which both angles swing at one common frequency \(\omega\), \(\theta_i = u_i\cos\omega t\). Substituting turns (10) into \(\Bigl[\tfrac gl\begin{pmatrix}2&0\\0&1\end{pmatrix} - \omega^2\begin{pmatrix}2&1\\1&1\end{pmatrix}\Bigr]\begin{pmatrix}u_1\\u_2\end{pmatrix} = 0\). Unless the determinant of that bracket vanishes, the only answer is \(u_1 = u_2 = 0\), the pendulum hanging still; setting it to zero gives \(2\,(g/l - \omega^2)^2 = \omega^4\), so

The first row of (10) then fixes the shape of each motion, \(u_2/u_1 = 2\,(g/l-\omega^2)/\omega^2\). For the slow mode this is \(+\sqrt2\): both arms swing together, the lower one about 1.4 times as far. For the fast mode it is \(-\sqrt2\): the arms swing against each other (Figure 2). These are the normal modes, and any small motion is a mixture of the two. Because \(\omega_+/\omega_-\) is irrational, the mixture never repeats exactly, but it stays perfectly regular.

06Chaos

At larger angles the small-angle picture breaks down, and the nonlinear sine and angular-speed terms matter. Some starts then show chaotic motion, although high energy alone does not guarantee it. The clearest sign is sensitive dependence on the starting state. Start two pendulums whose upper angles differ by \(0.001\ \mathrm{rad}\) and follow the distance \(d(t)\) between their lower masses. Figure 3 shows it for two starts, computed with the same equations and the same time step as the page.

Why do the two part? In small swings the equations are linear, so the difference between two motions is itself a small swing that keeps its size: that is the flat low-energy curve. In large swings the lower arm is carried up again and again to near the top, where it balances like a pencil on its tip, and a hair’s difference decides whether it goes over or falls back. Each such moment magnifies the gap, so it grows in bursts, most strongly while the lower rod is near the top.

A straight line on a logarithmic scale means \(d(t) \approx d_0\,e^{\lambda t}\). For the high-energy start of Figure 3 the curve wobbles, but a straight line laid through it, up to the moment the distance reaches \(0.1\ \mathrm{m}\) after \(4.1\ \mathrm{s}\), gives \(\lambda \approx 1.4\) per second: the distance doubles roughly every half second. The proper calculation follows a tiny change of the whole state, both angles and both angular velocities, shrinks it back whenever it grows large, and averages the logarithm of its growth over a long time. The limit of that average is the largest Lyapunov exponent,

where \(d\) is now the distance in the full state. Done in full for this start, with the page’s equations and step, shrinking the change back every \(0.1\ \mathrm{s}\) over \(20\,000\ \mathrm{s}\) of motion, it gives \(\lambda \approx 1.45\) per second; the 1.4 read off Figure 3, from one stretch of time and two of the four coordinates, was not far off. Its meaning is practical. Every extra second we want to predict costs a factor \(e^{1.45}\approx 4.3\) in the precision of the starting angles: to see ten seconds further, we would need to know the start about \(e^{14.5}\approx 2\times10^{6}\) times more precisely. Note what chaos is not: it is not noise, and it is not a failure of energy conservation. The energy stays fixed while the trajectory wanders over much of what that energy allows. Nor does a small energy drift guarantee more than itself: countless different paths share the same energy, so a drift of a few parts in \(10^{8}\) shows that the integrator is sound, not that the computed positions are the true ones beyond the predictability horizon.

07Solving it step by step

Equations (9) give the accelerations for any state, but no formula gives the angles at a later time. The computer therefore walks forward in small steps. Write the state as \(s = (\theta_1,\theta_2,\omega_1,\omega_2)\) and the equations as \(\mathrm{d}s/\mathrm{d}t = F(s) = (\omega_1,\omega_2,\alpha_1,\alpha_2)\). The fourth-order Runge–Kutta method looks at the slope four times per step, at the start, twice in the middle and once at the end, and takes a weighted average:

The error made in one step shrinks like \(h^5\) and the error accumulated over a fixed time like \(h^4\): halving the step makes the result sixteen times more accurate. The page uses \(h = 1/600\ \mathrm{s}\), ten steps per frame at sixty frames a second.

• Each frame advances by the real time elapsed multiplied by the Time scale, in steps of \(h\), with a 600-step guard. The clock caps a frame at 0.05 s and the slider caps the Time scale at 1.5, so even including the carried remainder no more than 46 steps can be requested; the 600-step guard cannot act.
• RK4 does not conserve energy exactly, which is why the Drift readout is there. It measures the signed difference \(E - E_0\) against the depth of the potential well, \(E_\text{scale} = (m_1+m_2)\,g\,l_1 + m_2\,g\,l_2\), rather than against \(E_0\), which can be close to zero. Over the first 20 s of the default start it stays within \(\pm 5\times10^{-8}\): its largest value, step by step, is \(+5.0\times10^{-8}\), at \(t \approx 8.7\) s; the readout, which samples once a frame, shows about \(+4.7\times10^{-8}\) there. It does not stay there. RK4 leaks energy, very slowly but steadily, so if you leave the default start running the drift creeps downward almost in a straight line, by about \(7\times10^{-8}\) a minute: after ten minutes it reads about \(-7\times10^{-7}\). Nor is the step equally good everywhere. Put the mass ratio at its top value of 10: each time the rods pass within a few degrees of a straight line, the heavy bob gives the light upper rod a sharp kick, with angular accelerations up to about \(1300\ \mathrm{rad/s^2}\) against \(75\ \mathrm{rad/s^2}\) for equal masses, and steps of \(h\) are too coarse to follow it. The drift then reaches \(-3.6\times10^{-4}\) within 20 s.
• Moving a mass, length, angle or angular-velocity slider restarts the motion from those values; the Time scale slider does not.

08From metres to pixels

The pivot \((x_0, y_0)\) sits at the centre of the figure, \(W\) pixels wide and \(H\) high, and the scale is chosen so that the fully stretched pendulum just fits: \(\sigma = (\min(W, H)/2 - 34)/(l_1 + l_2)\) pixels per metre. Screen \(y\) points down, so the minus signs of (1) turn into plus signs: bob 1 is drawn at \((x_0 + \sigma l_1\sin\theta_1,\ y_0 + \sigma l_1\cos\theta_1)\), and bob 2 a further \((\sigma l_2\sin\theta_2,\ \sigma l_2\cos\theta_2)\) beyond it. The protractor marks \(\theta\) every \(30^\circ\) from the downward vertical, on a faint circle of radius \(\sigma(l_1 + l_2)\), the farthest bob 2 can reach; the dashed circle inside it, of radius \(\sigma l_1\), is the track of bob 1.

The burgundy trail of the lower mass is drawn on its own transparent layer, one segment per RK4 step, and dimmed every frame by a factor \((1-0.008)^{60\,\Delta t}\), where \(\Delta t\) is the real time since the last frame. Its brightness would therefore halve every \(\ln 2/(60\,|\ln 0.992|)\approx 1.4\ \mathrm{s}\) of real time. The canvas, though, keeps opacity in whole steps of 1/255, and the small fading steps round away: the trail falls to half in about 1.6 s and then stops at about a quarter of its opacity, leaving a faint haze wherever bob 2 has been until Clear trail or a restart. Bob 1 has radius 8 screen pixels. Bob 2 has radius \(8(m_2/m_1)^{1/3}\) screen pixels, as for balls of the same material, clamped to 4–17 screen pixels.

09What the model leaves out

• Real rods have mass and a moment of inertia. A double pendulum made of two bars obeys equations of the same kind, with different coefficients.
• There is no friction. A real pendulum slowly loses energy, and its chaos eventually calms into small, regular swings.
• The masses may pass through each other and through the pivot; real rods would collide.
• After some fifteen to twenty Lyapunov times (\(1/\lambda \approx 0.7\ \mathrm{s}\) each; 11 to 14 s for the default start) the computed path is no longer the exact path from that start; it is a faithful example of how the pendulum can move. Averages over many swings are usually robust, but they deserve a check too: halve the step, lengthen the run, and see that they stay.

9.1Further reading

1. J. R. Taylor, Classical Mechanics, University Science Books, 2005: chapters 7 (Lagrange’s equations), 11 (coupled oscillators and normal modes) and 12 (nonlinear mechanics and chaos).
2. L. D. Landau and E. M. Lifshitz, Mechanics, 3rd edition, Pergamon, 1976: §5, problem 1, the plane double pendulum.
3. S. H. Strogatz, Nonlinear Dynamics and Chaos, 2nd edition, Westview Press, 2015.